function out = utilGains(PMI_RAW_T,RETURNS_T,SPREADS_T,RF_T,wLb,wUb,varMethod,preWindow,RiskAvCoef)

%Utility gains

numk = size(RETURNS_T,2);
T = size(RETURNS_T,1);

weights_Av = NaN(T,numk);
weights_NOSWITCH = NaN(T,numk);
weights_UNREST = NaN(T,numk);
PFRET_Av = NaN(T,numk);
PFRET_NOSWITCH = NaN(T,numk);
PFRET_UNREST = NaN(T,numk);
PFVar = NaN(T,numk);

for t = preWindow:T-1
        
    REC_1_to_t = (PMI_RAW_T(1:t,1) < 44.5)*1;
   
    RETURNS_2_to_t = RETURNS_T(2:t,:);
    SPREADS_1_to_t_1 = SPREADS_T(1:t-1,:);
    REC_1_to_t_1 = REC_1_to_t(1:end-1);   
    
    RETURNS_t1 = RETURNS_T(t+1,:);
    SPREADS_t = SPREADS_T(t,:);
    REC_t = REC_1_to_t(t);
    RF_t1 =  RF_T(t+1,:);
   
    
    RET_Av = mean(RETURNS_2_to_t);
    RETHat = NaN(1,numk);
    RETHat_Switch = NaN(1,numk);
    for j = 1:numk %loop over maturities

         %Newey-West std errors
         bv_j  = robustOLS(RETURNS_2_to_t(:,j),[ones(t-1,1) SPREADS_1_to_t_1(:,j)],0,0);
         bvSwitch_j = robustOLS(RETURNS_2_to_t(:,j),[1-REC_1_to_t_1 SPREADS_1_to_t_1(:,j).*(1-REC_1_to_t_1) ...
                                                     REC_1_to_t_1 SPREADS_1_to_t_1(:,j).*REC_1_to_t_1],0,0);

         RETHat(1,j) = [1 SPREADS_t(j)]*bv_j;
         RETHat_Switch(1,j) = [1-REC_t SPREADS_t(j)*(1-REC_t) REC_t SPREADS_t(j)*REC_t]*bvSwitch_j;
         
    end
       
    switch varMethod
        case 'Full'
            varHat_Full = var(RETURNS_2_to_t);
            varHat = varHat_Full;
        case '10y'
            varHat_10y = var(RETURNS_2_to_t(end-120+1:end,:));
            varHat = varHat_10y;            
        case '5y'
            varHat_5y = var(RETURNS_2_to_t(end-60+1:end,:));
            varHat = varHat_5y;         
        case 'EWMA'
            lambda = 0.97; %RiskMetrics monthly recommended value
            e = RETURNS_2_to_t - ones(t-1,1)*mean(RETURNS_2_to_t);
            vec_lambda = lambda.^(t-2:-1:0)';
            e_lambda = e.*(sqrt(vec_lambda)*ones(1,numk));            
            varHat_EWMA = diag((1-lambda)*(e_lambda'*e_lambda))';
            varHat = varHat_EWMA;
    end

    PFVar(t,:) = varHat;
    
    weights_Av(t,:) = 1/RiskAvCoef*RET_Av./varHat;
    weights_NOSWITCH(t,:) = 1/RiskAvCoef*RETHat./varHat;
    weights_UNREST(t,:) = 1/RiskAvCoef*RETHat_Switch./varHat;
        
    weights_Av(t,:) = max(wLb,min(weights_Av(t,:),wUb));
    weights_NOSWITCH(t,:) = max(wLb,min(weights_NOSWITCH(t,:),wUb));
    weights_UNREST(t,:) = max(wLb,min(weights_UNREST(t,:),wUb));
        
    PFRET_Av(t+1,:) = ones(1,2)*([weights_Av(t,:); ones(1,numk)].*[RETURNS_t1; RF_t1]);
    PFRET_NOSWITCH(t+1,:) = ones(1,2)*([weights_NOSWITCH(t,:); ones(1,numk)].*[RETURNS_t1; RF_t1]);
    PFRET_UNREST(t+1,:) = ones(1,2)*([weights_UNREST(t,:); ones(1,numk)].*[RETURNS_t1; RF_t1]);
   

end

util_Av = PFRET_Av(preWindow+1:T,:)-RiskAvCoef/2*(PFRET_Av(preWindow+1:T,:)-ones(T-preWindow,1)*mean(PFRET_Av(preWindow+1:T,:))).^2;
util_NOSWITCH = PFRET_NOSWITCH(preWindow+1:T,:)-RiskAvCoef/2*(PFRET_NOSWITCH(preWindow+1:T,:)-ones(T-preWindow,1)*mean(PFRET_NOSWITCH(preWindow+1:T,:))).^2;
util_UNREST = PFRET_UNREST(preWindow+1:T,:)-RiskAvCoef/2*(PFRET_UNREST(preWindow+1:T,:)-ones(T-preWindow,1)*mean(PFRET_UNREST(preWindow+1:T,:))).^2;

utilGains = [mean(util_NOSWITCH-util_Av)' mean(util_UNREST-util_Av)']*1200;

dmUtilGainsTstat = NaN(numk,2);
dmUtilGainsPval = NaN(numk,2);
for j = 1:numk
    
[dm_NOSWITCH,sedm_NOSWITCH] = robustOLS(util_NOSWITCH(:,j)-util_Av(:,j),ones(T-preWindow,1),2,1);
[dm_UNREST,sedm_UNREST] = robustOLS(util_UNREST(:,j)-util_Av(:,j),ones(T-preWindow,1),2,1);

dmUtilGainsTstat(j,:)= [dm_NOSWITCH/sedm_NOSWITCH dm_UNREST/sedm_UNREST];
dmUtilGainsPval(j,:)= 1-normcdf(dmUtilGainsTstat(j,:));

end

out.utilGains = utilGains;
out.dmUtilGainsTstat = dmUtilGainsTstat;

end
